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The  Generation  of  Three-Dimensional  Body-Fitted 
Coordinate  Systems  for  Viscous  Flow  Problems 

by  Z.U.A.  Warsi* 


Abstract 


'  The  main  aim  of  this  research  has  been  to  develop  and  implement 
a  technique  for  the  generation  of  spatial  coordinates  in  3D  regions 
enclosed  by  arbitrary  smooth  surfaces  for  ultimate  use  in  the 
numerical  solution  of  the  Navier-Stokes  equations.  In  this  regard,  a 
mathematical  model  based  on  a  set  of  elliptic  PDF’s-  has  been 
developed,  which  has  been  used  to  generate  smooth  coordinates  in  the 
region  formed  by  arbitrary  inner  and  outer  surfaces  of  known  shapes, 
around  multibodies,  particularly  around  a  wing-body  combination. 
These  equations  have  also  been  used  to  generate  surface  coordinates 
in  arbitrary  surfaces  and  are  also  capable  of  coordinate 


redistribution  in  any  desired  manner  both  in  3D  regions  and  in  2D 


surface  regions. 


"Principal  Investigator. 


Introduction 


7 The  problem  of  numerical  grid  generation  in  surfaces  and,  in 

(  -  r  ■>  *  -  w  i  ■' 

--'7 

three-dimensional  configurations  through  elliptic  PDE’s^has  been 
pursued  under  this  grant  , period.-  Various  reports  and  journal 
publications  produced  under  |the  present  grant  have  explained  the 
mathematical  model  in  detail.  \  (Refer  to  a  list  of  publications 
produced  in  the  grant  period). > The  developed  mathematical  model  has 
been  programed  on  CRAY-1  and  has  been  tested  for  single  and  two-body 
configurations  enclosed  in  a  single  outer  boundary  prefer  to  Appendix 

C)  and  for  generation  of  coordinates  in  a  single  surface, (refer  to 

! 

Appendix  B).  The  main  elements  of  the  mathematical  development  are 
contained  in  Appendix  A. 
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A  NOTE  ON  THE  MATHEMATICAL  FORMULATION  OF  THE  PROBLEM 
OF  NUMERICAL  COORDINATE  GENERATION* 

By 

Z  U  A  WARSI 

*  Mississippi  Stair  University 

Abstract  A  set  of  second  order  partial  differential  equations  for  the  generation  of 
three-dimensional  grids  around  and  between  arbitrary  shaped  bodies  has  been  proposed. 
Thtve  equations  basically  depend  on  the  Gauss  equations  for  a  surface  and  have  been 
structured  in  such  a  way  that  an  automatic  connection  is  established  between  the 
succeedinggenerated  surfaces. 

The  vanishing  of  the  Riemann  curvature  tensor  has  been  u>ed  to  isolate  those  funda¬ 
mental  equations  which  every  coordinate  system  in  either  two-  or  ifuee-dimer.>:onaI 
Euclidean  space  must  satisfy. 

1.  Introduction.  The  problem  of  generating  spatial  coordinates  by  nuivciicil  —  tUvJs  is 
a  problem  of  much  interest  in  practically  all  branches  of  engineering  science  and  phsics. 
At  present  a  number  of  techniques  are  under  active  development  for  the  generation  of 
two-  and  three-dimensional  coordinates  around  and  between  bodies  of  arbitrary  shapes. 
Among  these  efforts  two  easily  discernable  groups  can  be  formed:  (i)  the  method?  based 
on  the  solution  of  Certain  PDE's,  preferably  of  the  elliptic  type,  and  (ii)  the  a^c^aic 
methods.  In  a  large  number  of  methods  in  the  first  group  a  set  of  inhomogeneous  Laplace 
equations  is  talers  as  the  basic  generating  system.  These  equations  aie  then  inver'ed  and 
solved  for  the  Cartesian  coordinates  in  the  transformed  plane.  BaMrd  on  this  line  of 
approach  which  started  with  the  work  of  Winslow  [IJ  some  very  practical  results, 
particularly  in  two  dimensions,  have  been  obtained  by  Thompson  et  al  (2j  and  others.1 
For  the  generation  of  plane  curvilinear  coordinates  some  authors  have  also  used  hyper¬ 
bolic  and  parabolic  systems  of  equations,  (3]  For  the  methods  in  the  second  group,  i.e., 
the  algebraic  methods,  refer  foI3j. 

In  this  paper  we  have  first  considered  the  foi-muJalion  of  a  3 D  grid  generation  scheme 
which  is  basically  dependent  on  the  Gauss’  equations  for  a  surface  In  this  scheme  a  series 

*  Rece'ved  Mi)  10,  1982  This  rena/ch  has  beec  supported  b>  the  Air  Forte  Office  of  Scientific  Re?e»-ch, 
under  giant  AFOSR  80-0 1 85  The  material  of  this  p.apet  i»  bared  oe  an  extended  p*pet  juboimed  at  • 
NASA  AFOSR  epon-ored  oonh, cnee  on  Dutneri,*!  grid  genera'-ios  held  at  Nashville.  Te-.ntssee  Apnl  13-16. 
19*2.  [6] 

1  Refer  lo  |3]  for  as  exteny  ve  biV.ogr«ph>  oo  the  jubjeci 
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of  surfaces  arc  generated  bated  on  the  given  data  of  arbitrary  shaped  inner  and  outer 
surfaces  The  method  has  also  been  structured  in  such  a  way  that  the  variation  of  the  third 
coordinate  from  one  generated  surface  to  the  next  is  fully  reflected  in  the  system  of 
generating  equations.  It  has  been  found  that  these  generating  system  of  equations  are 
simple  to  implement  numerically.  In  particular,  the  solution  of  the  proposed  equations 
tends  to  the  solution  of  the  Laplace  equations  in  the  transformed  plane  in  case  the  surface 
becomes  a  Cartesian  plane.  An  evact  solution  of  these  equations  for  the  case  of  three-di¬ 
mensional  curvilinear  coordinates  between  a  prolate  ellipsoid  and  a  sphere  has  been 
obtained. 

In  a  plane,  or  a  surface,  or  a  3 D  space  there  are  endless  possibilities  of  introducing 
either  orthogonal  hr  non-orthogonal  coordinates.  This  realization  imparts  a  sense  of 
arbitrariness  to  the  choice  of  the  method  to  be  used  for  coordinate  generation.  If  it  is  a 
priori  decided  that  the  method  should  be  based  on  solving  partial  differential  equations 
then  the  arbitrariness  in  the  selection  of  the  appropriate  equations  for  the  generation  of 
coordinates  becomes  a  prob.cm  to  be  resolved  In  Sec  3  of  this  paper  it  has  beer,  shown 
that  despite  this  arbitrariness  it  is  possible  to  uncover  certain  sets  of  equations  which  must 
invariably  be  satisfied  no  matter  which  equations  or  methods  have  been  used  to  generate 
the  coordinates.  For  a  detailed  discussion  of  the  methods  discussed  here  and  on  some 
numerical  results  refer  to  Warsi  [4-7], 

2.  Gem  rating  svsft-m  ba-ed  on  the  Gauss*  equations.  In  the  process  of  formulation  of  a 
3 D  coordinate  generation  problem  it  i>  helpful  to  imagine  the  coordinates  of  a  point  in 
space  as  the  intersection  of  three  distinct  surfaces  on  each  of  which  one  coordinate  is  held 
fixed.  Using  the  convention  of  a  right-handed  curvilinear  coordinate  system  jr\  x:.  ;r  or 
(.  ij.  f  (refer  to  Fig.  1)  we  introduce  the  coordinates  in  the  surface  x*  -  const  through  the 
follow  ing  scheme. 

l:(*J,xJ);  r  =  2:  (*>.*');  r  =  3:  (x\  xJ) 


7 


Fig  I  RjJ,V  handed  coord  natc  arranfctmfm  and  basis  vectors 
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Thus  the  unit  norma!  vector  on  the  surface  x'  —  const  is 

rf'>=  (r<tXr,)/|r.  *r>|  (l) 

where 

r=  1:«  =  2./?  =  3;  w  =  2:  a  =  3,  fi  =  l  ;  p  =  3:  a  =1,0  =  2.  (2) 

From  elementary  differential  geometry  [8]  we  have  the  result  that  the  rectangular 
Cartesian  coordinates  r  =  (x,  y ,  t)  or  (x,,  x2,  x3)  of  any  point  on  every  surface  em¬ 
bedded  in  an  Euclidean  Ey  must  satisfy  the  equations  of  Gauss.  The  Gauss  equations  for  a 
surface  x'  =  const,  are  given  by 

=  Tir,  +  (3) 

where  all  the  Greek  indices  except  p  can  assume  only  two  values.  The  values  of  a,  0  and 
the  range  of  5  with  v  follow  the  scheme  given  in  (2)  In  Fq.  (3), 

_  3r  _  32r 

~  3x*  *  ~  3x'3x* ' 

T%0  are  the  surface  Christoffel  symbols  of  the  second  kind,2  i.e., 

r:fl  =  gmt[aP.o],  (4a) 

f  »  „1  -  1  (  a  9gf 

and  b ^  are  the  coefficients  of  the  second  fundamental  form.  Since  on  the  surface 
x'  =  const,  the  vector  is  orthogonal  to  the  surface  vectors  ra, 

^  =  rfr)r„.  (5) 

For  the  purpose  of  a  clear  notation  we  denote  the  space  Christoffel  symbols  as 

r,j  =  gtVj*k  ]-  (M 


where 


Using  (6a),  we  have 


r  =  —S-l—  =  r'r 

,J  Sr 'Sr/  ,J  " 


"  3x'3x;  "  ' 

where  all  the  Latin  indices  assume  three  values. 

To  Ex  ideas,  we  envisage  a  surface  which  is  formed  of  the  surface  coordinates  ((,  ij)  and 
on  which  {  -  const.  Dropping  the  index  v  in  Fq.  (3),  the  three  equations  for  the  second 
derivatives  of  the  Cartesian  coordinates  are 

=  Tf, r,  +  6„n,  =  Tfar,  4  bl2n,  r„  =  T'  r,  +  b*n.  (8a,b.c) 

In  Eqs.  (8)  a  is  orthogonal  to  both  and  rf  and  the  coefficients  of  the  first  fundamental 


*  Refer  tc  Appendix  A  for  t  eoIJectioo  of  other  formulae 
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fonn  gn,  g jj,  gn  are  assumed  lo  have  been  evaluated  at  l  =  const.  Obviously 

m 

'  *11  =  *\  +  >/  +  «*.  *12  =  *f*,  +  -Vf-V,  +  V  *22  =  -*,5  +  (9) 


■  ’  Multiplying  Eq.  (8a)  by  gjj,  Eq.  (8b)  by  -2 gl2,  Eq.  (8c)  by  g,,  and  adding  the  three 
equations,  we  get 

*+[(M)'<+(M)rjGj  =  n*,  (10) 

where 

^  =  *223{f  —  ^g,^,  +  *n3*v  (11*) 

A2  =  [3t{(*223l  ~  *!29,)//^3  )  +  3,{(*||3,  ”  *il3f)/v^3  (Ht>) 

^3  =  *11*22  (*12)  •  (He) 

n  =  (*,  y.Z), 

x  =  (y^i  ~  JV<)//gT,  y  =  (-V*  -  xf *,  Wcj.  2  =  (x(y,  -  X,yf  )//c7, 

Old) 

^  =  (*22^11  _  2g|2^)2  +  ^11^22)  ~  ^3(^1  "1"  ^2)’  (lit) 

where  *,  +  *2  is  twice  the  mean  curvature  of  the  surface. 


The  operator  A (V  or  simply  A2  defined  in  (lib)  is  the  second  order  differential  ope--ior 
of  BeJtranu  /8J.  for  the  surface  £  =  const  For  any  surface  x’  -  cons;  and  following  the 
scheme  (2),  we  have 

W =[a.{0toa.  -  *.,a,W£)  +  %{(*.A  -  *.A)//^)]/v'^,  02) 

where  G,  are  defined  in  Eq.  (A. 9).  It  is  easy  lo  show  by  using  the  surface  Christoffel 
symbols  that 

Aj{  =  (2g|jT|'2  -  g22T,',  —  *i|T^)/Gj,  (l?a) 

A2rj  =  (2g|2T,22  -  g22T|2,  -  *hT^j)/Gj,  (I3b) 

where  the  metric  coefficients  g0^  are  those  as  defined  in  (9).  It  is  interesting  to  see  that 
when  (he  Lapfacian  operator  v2  for  a  two-dimensional  Cartesian  space  is  applied  lo  the 
curvilinear  coordinates  (£,  rj)  in  an  Euclidean  plane,  the  resulting  expressions  are  exactly 
of  the  same  form  as  (1 3a,  b),  that  is  (refer  to  Eq.  (A.  13)), 

(2gl2r,'2  -g^r,1,  -*„ri)/(/)J,  (»<•) 

vN  =  (2g|2rjj  -  *22r2,  —  gnTy)/  (y) ,  ( Mb) 

where  now 

*n  =  xl  +  •  *i2  =  +  Ft>,.  *22  =  *1  +  >,*.  ^  =  *!>',  “  *,>r 

Though  the  right-hand  side  term  A  defined  in  (lie)  has  the  necessa-y  extrinsic  effects, 
nevertheless  we  must  have  an  explicit  dependence  of  r  -  (x,  y,  r)  on  the  third  coordinate 
{.  Thus  using  Eq.  (A  ll)  we  have 
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Ttl  =  I'i*  1  rg  +  I'm1-,  +  I'  |3 1  rf , 

(15a) 

r(^  —  I'l2r<  "*■  r*  7"  r t32  fj , 

(15b) 

r«  =  +  Yi2rn  +  , 

(15c) 

and  we  evaluate  each  of  these  derivatives  at  {  =  const.  TaJu'ng  the  dot  product  of  Eqs  (15) 
with  n  and  comparing  with  Eqs.  (5),‘  we  get 

h,,  —  Xrf|.  =  Xr,J2,  bj2  —  Xr22,  (16a) 

where 


X  =  d  •  r{  =  Xx(  +  Yy(  +  Zz( .  (16b) 

Thus,  the  expression  (1  le)  for  R  is  replaced  by 

^  =  +  (17) 

2.1  Fundamental  generating  system  of  equations.  We  now  impose  the  follow  ing  differen¬ 
tial  constraints  on  the  coordinates  £  and  17: 

A2C  =  0.  A2tj  =  0,  (18) 

and  take  them  as  the  basic  generating  equations  for  the  coordinates  in  a  surface  A 
comparison  of  Eqs.  (13)  and  (14)  has  alreadv  shewn  that  is  not  a  2D  Laplace  operator 
except  when  the  surface  degenerates  into  a  plane  having  no  dependence  on  the  .* -coordi¬ 
nate. 

It  is  a  well-known  result  in  differential  geometry  that  the  i'-cthermic  coordinates  in  a 
surface  satisfy  Eqs.  (18)  identically.  The  isothetmic  coordinates  t  and  tj  are  those 
orthogonal  coordinates  in  a  surface  winch  yield  gu  =  gn.  The  situation  here  is  parallel  to 
the  choice  of  the  Laplace  equations  =  0,  V3n  =  0  for  the  generation  of  plane 
curvilinear  coordinates  (2],  which  are  also  satisfied  by  the  confoimal  coordinates -in  a 
plane.  The  important  point  to  note  here  is  that  the  satisfaction  of  the  Laplace  equations  is 
a  necessary  but  not  a  sufficient  condition  for  the  existence  of  conformal  coordinates. 
Similarly,  the  satisfaction  of  equations  (18)  is  a  vc-.e'-xuy  but  not  a  sufficient  condition  for 
the  existence  of  isothermic  coordinates.  It  would,  therefore,  be  more  meaningful  if  we 
interpret  Eqs  (18)  as  providing  a  set  of  differential  constraints'  on  the  metric  coefficients 
gn,  S 12.  a™3  iu  defined  in  (9). 

Having  chosen  Eqs.  (18)  as  the  generating  system,  the  equations  for  the  determination 
of  the  Cartesian  coordinates,  viz.,  Eq.  (10)  becomes 

Cr^nR,  (19) 


where  £,  n  and  R  have  been  defined  in  (11a),  (1  Id),  and  (17)  respectively.  The  three  scalar 
equations  in  expanded  form  are  now 

~  +  *11*,,  =  XR’ 


~  Zguhf,  +  fnv„  =  YR, 


1 A  m unites u ti on  of  the  rr.xnv  possibilities  for  introducing  coorAicxtes  in  •  given 


(20a) 

(20b) 
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g22*i(-  2*i2*<,  +  fn*„  =  ZR.  (20c) 

For  a  plane  z  =  const.,  R  =  0  and  the  Eqs.  (20)  are  the  inversions  of  the  Laplace 
equations  in  the  ^rj-plane. 

It  can  be  shown  that  Eqs.  (20)  can  be  combined  to  obtain  the  equations  of  a  surface 
z  —  z(x,  y)  in  the  well  known  form, 

azxx  -  2$zxy  +  yzyy  =  2HM,  (21) 

where 

2H  =  k,  +  k3  =  R/Gy,  M=\+pi  +  q7,  p  =  zx,q  =  zy, 

a  =  (1  +  ?)/&,  P  =  pq//M,  y  =  (l  +  p7)/fi . 

Using  the  following  definitions  and  identities 

=  («.2)2-  X=-p/{M,  Y  =  -q/ ,  Z=\/,fM, 

Hl(x,x)  =  (\  -X7)G>,  A,(*,y)  =  -XYGit  A,(>.  y)  =  (1  -  Y7)G3, 

where 

*\(a>b)  =  +  o,*t)  + 

calcu'ating  zf{,  z£l),  from  z£,  z,  and  substituting  these  expressions  in  (20c)  while  using 
Eqs.(20a,b),  we  get  Eq  (21). 

We  new  use  the  result  that  if  (£.  »j)  is  a  permissible  system  of  coord- .nates  in  a  surface 
then  so  is  (l,  rj),  where  |  =  £(£,  ij ),  rj~v((,v),  provided  that  the  Jacobean  of  the 
transformation  does  not  vanish.  It  is  a  straight  forward  matter  to  shew  that  on  coordinate 
transformation,  Eqs  (20)  become 


fjc  =  XR,  Ly  =  YR,  £z  =  ZR, 

(22a.  b.c) 

*  =  82 2*i\  ~  2EAa  +  fn^,-  +  +  Qdj- 

(23a) 

P  -  iu^u  ~  ^8)2^12  +  81 1^22' 

(23b) 

Q  ~  822? u  ~  ^8\2^\2  +  8\\Pn< 

(23c) 

,y  _  dx •  Bx*  d7r 

49  a***  a*9  a x*dx*' 

(23d) 

and  X,  Y,  Z,  and  R  are  exactly  the  same  expressions  as  given  in  (1  Id)  and  (17),  in  the  xm 
coordinate  system  It  is  preferable  to  solve  Eqs.  (22)  with  P£t  as  arbitral!)  prescribed 
functions  of  the  coordinates  This  aspect  of  the  formulation  therefore  provides  a  capability 
to  redistribute  the  coordinate  systems  in  the  surface  in  any  desired  manner.4 

2  2  Example  of  a  solution  algorithm  The  discussion  that  follows  pertains  to  the  case 
when  it  is  desired  to  generate  the  3 D  curvilinear  coordinates  between  two  arbitrary  shaped 
smooth  surfaces  Let  the  surface  coordinates  (f,  i)  of  the  inner  body  ij  =  ij t  and  of  the 

4  For  *  limiting  form  of  Eqs  (22)  refer  to  Appendix  B. 
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outer  bod>  ij  =  »?*  be  the  same  coordinates.  Thus 

V  =  Vs  r  =  rt(t,{);  =  »!.:  r  = 

are  known  functions  (either  analytically  or  numerically)  and  hence  the  needed  partial 
derivatives  with  respect  to  t  and  £  are  directly  available  at  the  surfaces. 

For  the  computation  of  rf  in  the  field  one  must  first  note  that  the  coordinate  £  may  not, 
in  genera],  satisfy  the  equation  A'f’f  '=  0.  Consequently,  must  satisfy  the  equation 

iP'r  +  G2(^)ts  =  G2(*P'  +  AfV2'-  (24) 

From  this  equation  we  devise  a  weighted  integral  formula5 

*/[/, (l)frtt),  +  fiU)( r„)  J  (25a) 

where 


*n  (dM  /g2  / 


(25b) 


and 

/i(»»*)=l.  /,(’»«)  =  0.  f2(vB)~  0.  /2(i?*)  =  l. 

rhere  is  no  difficulty  in  the  numerical  evaluation  of  (25a)  in  an  iterative  cycle  after  the 
weighting  functions /,  and  f2  have  been  prescribed  a  priori. 

Referring  to  Fig  (2a).  we  now  solve  Eqs.  (20)  or  (22)  for  each  l  ~  const.,  by  prescribing 
the  values  of  x,  y,  z  on  the  curves  C,  and  C2  which  respectively  represent  the  curves  on  B 
and  oc.  In  Fig.  (2b)  C}  and  C4  are  the  cut  lines  on  which  periodic  conditions  are  to  be 
imposed. 

2.3  An  exact  solution  of  the  proposed  equations.  The  following  example  demonstrates  that 
the  proposed  set  of  generating  equations  (22)  are  consistent  and  provide  nontrivial 
solutions. 

We  consider  the  case  of  coordinate  generation  between  an  inner  body  ij  =  *iB  which  is  a 
prolate  ellipsoid  and  an  outer  boundary  jj  =  r?x  which  is  a  sphere  The  coordinates  which 
vary  on  these  two  surfaces  are  (  and  A  curve  C,  on  the  inner  surface  for  S  ~  const,  is 

x  ~  t  cosh  ijB  cos  f ,  y  =  Tsinhrjisini’cos(,  z  -  rsinhij^sinfsin  £.  (26a) 

where  t  and  ijB  are  the  parameters  of  the  ellipsoid  Similarly  tv  r  curve  C2  corresponding  to 
the  same  {  =  const  on  the  outer  surface  is 

x  -  exp(7iot)cosf,  >■  -  exp(n*)rinfcos£,  z  -  exp(Tja  )sin  £  sin£.  (26b) 


JThf  discussion  jivro  b err  is  directed  to  the  stuabon  of  Fig  (2a)  For  other  situations.  eg  .  simptv-con- 
ceded  domains  or  multi  bod)  problems  the  method  of  calculating  mus'  aJ*ays  be  devised  separaiel)  Note 
also  «hat  Eq  (24)  reflects  the  condition  ifj'(  -  0 
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F)C.  2.  (»)  Topology  of  ihe  gi'fn  surfaces 

(b)  Surface?  to  bt  generated 


In  order  to  generate  surfaces  on  which  (  and  ij  are  the  coordinates  and  in  which  the 
ij-coordinate  can  nonuniformly  be  distributed  (c*  -  traction  or  expansion  in  the  tj  coordi¬ 
nates),  we  assume 

£  =  £(£).  V  =  v(ii)  +  Va>  (27) 

where  £  =  0  corresponding  to  £  ~  0  and  ij  =  rjB  corresponding  to  ij  =  i?,.  Thus  ((0)  =  0, 
T)(rjB)  —  0.  Under  the  transformation  (27),  the  only  nonzero  components  of  PJe  are 
and  Pj2.  Writing 


we  have 

Based  on  the  forms  of  the  boundary  conditions  (26a, b)  we  assume  the  following  forms 
for  x,  y ,  z  for  each  f  =  const.: 

x=f(ij) cos!,  y  =  $(Tj)sinfcos{,  i  ~  )sin£sinf.  (29) 

The  boundary  conditions  are 

/(la)  =  Tcoshri,.  /(vm)  =  t'Pfv*). 

<»(^*)  =  Tsinhn,,  4>(rjm)  =  expfij,).  (30) 

Using  the  expressions  of  (29)  we  calculate  the  various  partial  derivatives,  metric  coeffi¬ 
cients,  and  all  other  data  as  needed  for  Eqs.  (22).  On  substitution  we  get  an  equation 
containing  sin*f  and  cosJf.  Equating  to  zero  the  coefficients  of  sinlf  and  cosJf,  we 
obtain 


•*1  "Si 

II 

««>-&• 

1 

02  1  dB 

X  d\' 

f»— (ii 

(28) 
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/"//'  =  e/e  +  (3i) 

=  *'/*  +  ♦'/♦.  (32) 

where  a  prime  denotes  differentiation  with  respect  to  n  On  direct  integrations  of  Eqs.  (31) 
and  (32)  under  the  boundary  conditions  (30),  we  get 

f(v)  =  ^exp(2?ij(ij))  +  C,  (33a) 

<f>(i  )  =  Dt\p(BT}(ij)),  (33b) 

where 

A  =  r[(exp(ijoe)  -  -rcoshnjJsinhijJ/Cexpfij,,.)  -  rsinhij*), 

B  -  0l»  ~  ln(.rsinhTj#  ))/(i|.  -  ij,), 

C  =  r[exp(T?x)(coshij,  -  sinh^)]/(exp(ijx)  -  rsinhij^), 

D  =  r  sinh  ^ . 


As  an  application,  we  take 


1(1)  =  oh 

where  a,  b,  and  k  are  constants.  Thus 

n(i)  =  t>(i ~  Vb)*'. 

(34) 

rfO-kv 

-  ’»*)(’? VeKa-i.) 

(35) 

By  taking  a  value  of  k  slightly  greater  than  one  ( k  -  1.05)  we  can  have  sufficient 
contraction  in  the  ^-coordinate  near  the  inner  surface.  For  the  chosen  problem  since  the 
dependence  on  f  is  simple,  we  find  that  the  generated  coordinates  between  a  prolate 
ellipsoid  and  a  sphere  are 

x  =  [A  ))  +  C]cosf,  y  —  Z?exp(2?rj( rj  ))sLn  fcosf. 

2  =  f>exp(Brj(^  ))sin  f  sin  {.  (36) 

This  example  also  shows  that  the  chosen  generating  system  of  equations  (20)  or  (22)  arc 
capable  of  providing  non-isothermic  coordinates  h-etween  smooth  surfaces. 


3.  Differential  equations  hived  on  the  Riemann  tensor.  In  any  given  space  there  are 
endless  possibilities  for  the  introduction  of  coordinate  curves  Each  chosen  set  of  curves 
determines  its  own  metric  components.  For  example,  in  a  Cartesian  plane  besides 
introducing  rectangular  Cartesian  coordinates  x,  y,  we  also  have  endless  possibilities  foT 
introducing  either  orthogonal  or  noncrthogonal  coordinate  curves.  However,  as  is  well 
known,  there  is  a  basic  differential  constraint  on  the  variations  of  g,* s  irrespective  of  the 
coordinate  system.  Since  the  curvature  of  an  Euclidean  two-dimensional  plane  is  identi¬ 
cally  zero,  the  basic  differential  constraint  on  the  g  ’ s  is 


(GyY'^R 


a 

H 


(37) 


where  ij  are  any  arbitrary  coordinate  curves  in  the  plane.  Thus  no  matter  which 
coordinate  system  is  introduced  in  a  plane,  the  corresponding  matrics  g(y  must  satisfy  Eq. 
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(37).  Equation  (37)  has  also  been  used  as  the  basic  generating  equation  for  the  generation 
of  orthogonal  coordinates  in  a  plane  [9]  In  general,  the  Riemann  curvature  tensor  #rjnf 
defined  as. 


-i(; 


jr, i _ j[  &JP  \ 

$xp  9x'dx"f 


9* '8*'  dxJdx'  9*'3x"/ 

+g"\[j"' *][*#'*]  -IjMII'".*!)  (38) 

defines  the  components  of  the  curvature  tensor  of  any  general  space.  If  the  space  is 
^-dimensional,  then  the  number  of  components  #,jnf  are  given  by  A2(A 3  —  1)/12.  Thus 
for  A  =  2  there  is  one  distinct  surviving  component  stated  in,  Eq.  (37).  However,  for 
A  =  3,  it  has  six  distinct  components  #l2t2-  #23231  #12131  #12321  # uu •  M  the 

3D-space  is  Euclidean,  then  its  curvature  is  zero,  so  that  the  six  equations 

#1212  ~  #1313  ~  0,  #2323  ~  ®» 

^1213  =  0.  ^UJ2  =  O'  ^1323  =  0  (39) 

determine  the  differential  constraints  for  the  six  metric  coefficients  g  in  any  coordinate 
system  introduced  in  an  Euclidean  space.  These  equations  in  the  expanded  form  have 
been  given  in  [5]  and  [6]. 

Equations  (39)  are  those  consistent  set  of  partial  differential  equations  which  must 
always  he  satisfied  by  the  metric  coefficients  gtJ.  In  the  3 D  case  Eqs.  (39)  are  six  equations 
in  six  unknowns,  and,  therefore,  they  form  a  closed  system  of  equations.  In  contrast,  for 
the  2D  case  there  is  only  one  equation  (Eq.  (37))  and  three  unknowns  gu,  gi2,  and 
therefore  some  constraints  are  needed  to  turn  Eq.  (37)  (such  as  orthogonality  [9])  into  a 
solvable  equation. 

To  obtain  the  Cartesian  coordinates  on  the  basis  of  the  available  gtJ’ s,  we  introduce  the 
unit  base  vectors  A,  as 

A,  =  a ,//g„  ,  no  sum  on  /.  (40) 

Let  the  components  of  A,  ak>ng  the  rectangular  Cartesian  axes  be  denoted  as  ut,  c,,  h,,  so 
that 


A,  =  (»,,  «>,,*;), 

where 

“1  =  x(//gu  •  ci  ®  V/*n  •  Hi  =  V^t  * 

«2  =  X*/\t!22  1  V2  ~  V,/V^22  .  *2  =  2  ' 

«3  =  .  *1  =  .  *j  =  *t/Jg33  ■  (4 1 ) 

If  the  components  u,,  v,,  w;  become  known  by  some  method  then  it  is  possible  to  evaluate 
the  Cartesian  coordinates  through  the  line  integrals 

r  =  +  ^2\fiudr\  +  ^3>^33^)‘  ^ 


(«) 
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Tbe  determination  of  r(f  wt  is  a  separate  problem  which  we  now  consider.  First  of  all 
using  (40)  in  Eq.  (A.l  1),  we  get  a  system  of  first  order  partial  differentia]  equations 


r2 

•j 


(43) 


where  there  is  no  sum  on  the  repeated  index  /'.  Equations  (43)  form  a  system  of  27  first 
order  PDE’s  in  nine  independent  variables  ut%  t>(,  w,.  This  system  of  equations  is 
oc-erdetermined  and  thus  its  solvability  should  depend  on  certain  compatibility  conditions. 
According  to  a  theorem  on  the  overdetermined  system  of  equations  [10],  if  the  compatibil¬ 
ity  conditions  bold  then  the  solution  of  Eqs.  (43)  exists  and  is  unique.  The  conditions 

B2Aydxmdx'  =  d2A,/dx'dxm  (44) 


for  all  values  of  »,  m,  and  j  are  the  compatibility  conditions.  To  prove  (44)  we  use  Eq. 
(A.  1 1),  which  on  cross  differentiation  yields 


dxmBx' 


dxJdxm 


(45) 


where  R'.,mj  is  the  Riemann-Christoffe?  curvature  tensor  and  is  related  to  the  Riemann’s 
tensor  Rljkl.  Evidently  in  our  present  case  R'.imj  =  0,  since  the  space  is  Euclidean. 
Inserting  (40)  in  (45)  wc  find  that  Eq.  (44)  axe  identically  satisfied. 

It  is  interesting  to  note  that  for  a  tw  o-dimensional  curvilinear  coordinate  system  there  is 
no  need  to  solve  the  system  of  equations  such  as  (43).  In  this  case  the  single  differential 
equation  with  G,  =  g 


implies  the  existence  of  a  single  function  a((,  tj)  such  that 

-Jg 

«,  =  -  -r2. 


&ii 


.  -  —  r> 

°-  —  1  12' 


* 


II 


Consequently 

n,  =  cosa,  u,  =  -sina,  u2  =  cos(a  -  9),  »2  = -sin (a  -  6), 

where  a  is  the  angle  made  by  the  tangent  to  the  coordinate  line  ij  =  const  in  a  clockwise 
sense  with  the  x  axis,  and 

cos  e  =  g[1//gug23 


is  known  The  angle  o  becomes  known  since  gtJ  are  known,  e  g.  [9J. 

3  1  Case  of  orthogonal  coord  notes.  For  orthogonal  coordinates  since  the  cosines  of  the 
angles  between  the  coordinate  curves  are  zero,  we  have 

f  12  =  in  =  tu  =  0. 


(46) 
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Consequently, 

[12, 3]  =  [13,2j  =  [23, 1]  =  0,  r,3j  =  r,J3  =  ri  =  O,  g  =  gnt22833- 

The  equations  for  the  metric  coefficients,  viz.  Eqs  (39)  under  the  constraints  of  ortho¬ 
gonality  (46)  simply  reduce  to  the  Lame’s  equations  They  can  concisely  be  written  as  six 
equations  by  dropping  the  summation  convention  in  the  form 

JL/±  JL/±  *5tW-L  M*  M*-0  (47a) 

3x'  ( 3x'  I  +  3x4  [  A*  3x*  /  hj  3x'  3x'  *  * 

- 1  M*  +  JL  M±  ,.7K v 


3aA,  _  3A,  3 hj  ]  3/r^ 

3x'3x*  “  ^  3^>  &r*  K  bTk 


where  (/,  /,  A:)  are  to  be  taien  in  the  cyclic  permutations  of  (1,  2, 3),  in  this  order,  and 

^1  =  tfSu  »  *2  =  2  ’  ^3  =  /?33  ■ 

To  obtain  the  differential  equations  for  the  Cartesian  coordinates  x,  =  x,  x2—  y , 
jt3  =  2,  we  first  proceed  from  Eq.  (A  13)  and  have 

A).  ^vJn  =  /-(/t,V*2).  r*vV=^(A(V^).  («) 

where 

4*  -  /r,h2h3,  V  2  =  3„  4  3„  +  3„. 

Proceeding  directly  from  Eq.  (A.  14)  and  using  Eqs.  (46)  and  (48),  the  equations  for  the 
Cartesian  coordinates  are 

Zxm  =  0,  m=  1,2,3,  (49) 

where 

*  =  A  /*A  A l  +  l/Mj  A)  +  A  (*i*i  A ) 

■  an  *,  nl  ail  h,  aW  art  *,  nr 

Note  that  the  operator  Z  and  the  Laplacian  V 2  are  rela'ed  as 

z  *  = 

where  ^  is  a  scalar. 

Equations  (47)  and  (49)  are  those  consistent  set  of  equations  which  every  orthogonal 
coordinate  system  must  satisfy. 

3.2  The  rase  of  isothermic  coordinates.  Isothermic  coordinates  in  a  surface  embedded  in  a 
3 D  Euclidean  space  Are  those  coordinates  in  which  the  metric  coefficients  g,,  and  g,}  in 
the  surface  jj  =  const  are  equal  That  is,  the  element  of  length  ds  on  ij  -  const  is  given  by 

(ds)\-^  =  gu[(dif  ^  (dl)% 

where  (,  {  are  chosen  to  be  the  surface  coordinates  Using  (46)  and  setting 

gfJ  =  g„  and  gjj  =  F(ij) 
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in  Eqs.  (39),  we  obtain  tbe  basic  equations  for  which  are 


It  can  easily  be  verified  that  the  only  solution  of  Eqs.  (50c,  d)  is 

f„  =[a  +  />(ij)]7(*.  f).  a  =  «>ns‘-  (51) 

Thus  from  (50b) 

F(t,)  =  (dP/duf.  (52) 

Substituting  (51)  and  (52)  in  Fq.  (50a),  the  differential  equation  for/(f.  f )  becomes 

A|]i/|+i(l^|+2/=0.  (53) 

ad/3fran/af/  ; 

In  Kreyszig  [11],  we  have  the  result  that  if  in  a  portion  of  a  surface  isothermic 
coordinates  can  be  introduced  then  that  portion  of  the  surface  can  conformally  be 
mapped  onto  a  plane.  Thus  in  effect  the  solution  of  Eq.  (53)  provides  that  mapping 
function  which  conformally  maps  a  surface  onto  a  plane  As  a  verification  of  the  above 
conclusion,  we  verify  that  tbe  function 

/  =  4e,f/(I  +  e21)1  (54) 

is  a  solution  of  Fq.  (53).  This  function  is  related  with  the  isothermic  coordinates  on  a 
sphere.  Using  the  parametric  equations  of  a  sphere 

X  =[fl  +  /,(rj)]cos0,  y  =[a  +  F(»j)]sin5sin*,  i  =[a  +  />(i?)j>in0cos* 

and  writing 

9 

£  =  *.  f  =  In  tan  j , 
where  0  <  $  <  2w  and  0  <  8  <  v,  we  obtain 

g»  =  !„=«(«  4  + 

Thus  the  equations 

x  =  (a  +  P)(  1 -e’f)/0 
y  =  2(o  +  P)e'sint/0  +  <»). 

*  =  2(a  +  /^efcos^/O  +  en) 


(55) 
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represent  a  sphere  of  radius  a  ■+  ^(ij)  in  terms  of  the  isothermic  coordinates  £,  f  in  the 
surface  Since  P(i i)  is  an  arbitrary  function  of  jj,  we  now  have  the  capability  of  prescribing 
a  suitable  function  P(i j)  to  achieve  any  sort  of  contraction  on  expansion  in  the  field.  It  is 
expected  that  the  representation  (55)  should  prove  useful  in  the  computational  problems 
associated  with  a  sphere. 

4.  Conclusions.  In  Sec.  2  of  this  paper  a  set  of  second  order  PDE's  have  been  obtained 
which  generate  a  series  of  surfaces  between  the  given  inner  and  outer  arbitrary  shaped 
bodies.  The  necessary  mathematical  apparatus  which  connects  one  generated  surface  with 
its  neighbor  along  with  the  curvature  properties  of  each  surface  has  been  incorporated  in 
the  right  hand  side*  terms  of  the  equations.  (Eqs  (20)  or  (22)).  By  changing  the  computa¬ 
tional  techniques  these  equations  can  also  be  used  to  generate  the  3 D  coordinates  when 
more  than  one  inner  bodies  are  present  in  the  field. 

In  Sec.  3,  based  on  some  basic  differential  geometric  concepts,  a  number  of  field 
equations  have  been  isolated  which  must  always  be  satisfied  by  any  coordinate  system  in 
an  Euclidean  space.  Efficient  numerical  methods  are  to  be  developed  to  solve  these 
quasilincar  equations  (Eqs.  (39))  on  a  digital  computer. 

Appendix  A.  In  this  appendix  we  collect  some  useful  formulae  which  have  been  used  in 
the  main  text. 

As  Doted  in  the  text,  a  general  curvilinear  coordinate  system  is  denoted  as  x\  i  =  1,2,3, 
or  as  £,  rj,  f,  while  a  rectangular  Cartesian  system  is  denoted  as  xm,  m  =  1,2.3  or  as 
x,  y,  z  Since  r  is  a  position  vector  in  an  Euclidean  space,  the  covariant  base  vectors  a,  arc 
given  by 

=  3r/3x\  (A  1) 

while  the  contravariant  base  vectors  a'  are  given  by 

g'  =  gradx’.  (A  .2) 

The  covariant  and  the  contravariant  metric  components  are  respectively  given  by 

*„  =  *,  *,.  S'J  =  *  (A-3) 

Both  metric  coefficients  are  related  through  the  equations 

g'Jg,k  =  *i’  (A4) 

where  8{  are  the  Kronecler  deltas.  Also 

g  =  det(g''),  gg  =  I.  (A5) 

Based  on  (A.4),  we  also  have 

=  g'i%t  (A  .6) 

=  »1'/h)/2£  (a?) 

where,  here  and  in  all  the  expressions  a  repeated  lower  and  upper  index  always  stands  for 
•  sum  over  the  range  of  index  values.  Also 
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*  =  dct(jv)  =  g33G3  +  gl3Gs  +  g2JG* 

=  822^*2  +  SliG 4  +  t22^f>  ~  Sll^i  "*■  §12^4  (A.8) 


where 

G\  =  iuBii  ~  (&2i)  •  =  Snia  ~  (fi>)  >  =  £ii#22  “  (£12 )  • 

^4  =  f  uf23  —  fuf»»  =  812823  ~  812822 ’  &t  =  812812"  8uBu-  (A-9) 

Note  that 


«"  =  C,/g.  «“  =  C7a/g,  gn  =  G3/g, 
g'^G^g,  gn  —  Gs/g,  gv  =  Gt/g. 
The  derivative  of  a  covariant  base  vector  is  given  by 

faydx'  =  d^r/dx’dx'  =  fy,. 

The  Laplacian  of  a  scalar  $  in  a  curvilinear  coordinate  system  is 


*  *  U*’3jc'  "8*7 


(A. 10) 
(All) 


(A  12) 


where  T'j  have  been  defined  in  Eqs.  (6). 

If  $  =  xm  is  any  curvilinear  coordinate  then  from  (A. 1 2) 

vV=-g"r,;.  (a  13) 

If  $  =  *„  is  any  rectangular  Cartesian  coordinate  theD  from  (A.  12) 

*''^ib  +  (*v,£  =  0-  (AM) 


Appendix  B.  In  niimerical  computations  it  is  desirable  to  solve  Eqs.  (22)  in  their  limiting 
forms  on  certain  special  lines  in  the  field.  Referring  to  Fig.  (2a),  let  the  x-axis  be  aligned 
to  pass  through  the  inner  body  from  two  of  its  points,  which,  when  extended  in  both 
directions  meets  the  outer  body  at  its  two  corresponding  points.  The  portions  of  the  lines 
between  the  inner  and  the  outer  bodies  form  the  right  and  the  left  segments.  On  each 
segment  y  -  1  =  0,  and  according  to  the  adopted  convention  f  =  0  and  (  =  *on  the  right 
and  the  left  segments  respectively  for  all  values  of  t.  With  this  choice  of  the  axes  only  Eq. 
(22a)  is  of  interest.  Talung  the  limit  of  Eq.  (22a)  as  £  —  0  or  f  -*  *r,  we  obtain 

*ii  +  (Bl) 

or 

t  — 

where  the  control  function  P £  bas  already  been  chosen  a  priori.  The  terms  on  the  right 
hand  side  of  Eq.  (B.I)  are  difficult  to  assess  tor  their  limiting  behaviors  However,  some 
guidance  can  be  obtained  from  the  exact  solution  discussed  in  Sec.  2.3.  This  approach 
suggests  that  in  any  case,  the  following  estimates  can  be  used. 
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•*’=/i(^)(Ih),/2A,-  *  =Mv)< 

?*,  =f)(n)(iu)'/3,  f22  =  *i.  forf-Oorir,  (B  2) 

where /i./j./j  are  functions  of  if  Using  the  estimates  (B  2)  in  Eq.  (B.l),  we  obtain 

+  An  )*,-  =  0,  (B.3) 

where 

The  scheme  now  is  to  solve  Eq.  (B.3)  by  prescribing  7\^)  arbitrarily  to  achieve  the 
desired  control  pf  points  on  the  segments.  Since  has  already  been  chosen  in  advance 
this  approach  produces  those  values  of  F{t})  (though  they  Deed  not  be  calculated)  which 
are  consistent  with  the  basic  equation,  viz.,  Eq.  (22a). 
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ABSTRACT 


In  this  paper  a  set  of  second-order  partial  differential 
equations  for  the  generation  of  coordinates  in  a  given  surface  have 
been  developed  and  then  solved  numerically  to  demonstrate  their 
versatality.  The  proposed  equations  are  not  some  arbitrarily  chosen 
equations  but  are  a  consequence  of  the  formulae  of  Gauss  and  thus 
carry  with  them  an  explicit  dependence  on  the  geometric  properties  of 
the  given  surface.  Furthermore,  these  equations  are  easy  to  solve 
and  require  only  the  specification  of  the  bounding  curves  to  provide 
the  Difichlet  boundary  conditions  for  their  solution.  Results  of 
coordinate  generation  both  in  the  simply  and  doubly-connected  regions 
on  some  known  surfaces,  with  the  option  of  coordinate  redistribution, 
have  been  presented.  Extension  of  this  technique  to  arbitrary 
surfaces  seems  to  be  straightforward. 
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REPRODUCED  AT  GOVERNME.v;r  expense  •. 

I.  INTRODUCTION 

The  problem  of  generating  spatial  coordinates  by  numerical 
methods  through  carefully  selected  mathematical  models  is  of  current 
interest  both  in  mechanics  and  physics.  A  review  of  various  methods 
of  coordinate  generation  in  both  two  and  three-dimensional  Euclidean 
spaces  is  available  in  [1],  and  reference  may  also  be  made  to  the 
proceedings  of  two  recent  conferences,  [2],  [3].  on  the  topic  of 
numerical  grid  generation. 

This  paper  is  exclusively  directed  to  the  problem  of 
generation  of  a  desired  coordinate  system  in  the  surface  of  a  given 
body  and  thus,  in  a  basic  sense,  it  is  an  effort  directed  at  the 
problem  of  grid  generation  in  a  two-dimensional  non-Euclidean  space. 
The  mathematical  model  selected  for  this  purpose  is  based  on  the 
formulae  of  Gauss  for  a  surface  and  nas  been  discussed  by  the  author 
in  earlier  publications,  [M]  -  [7].  The  resulting  equations  are 

three  coupled  quasilinear  elliptic  partial  differential  equations 
with  the  Cartesian  coordinates  as  the  independent  variables.  These 
equations  are  nonhomogeneous  with  the  righthand  sides  depending  both 
on  the  components  of  the  normal  and  the  mean  curvature  of  the 
surface.  These  equations  therefore  reflect  the  geometrical  aspect  of 
the  surface  in  an  explicit  manner. 

The  proposed  equations  have  been  used  to  generate 
three-dimens ional  coordinates  between  two  given  surfaces  by  using  an 
extrinsic  form  of  the  mean  curvature,  [7],  [8],  [9].  However,  if  the 
purpose  is  to  generate  the  coordinates  only  in  a  given  surface  then 
the  intrinsic  form  of  the  mean  curvature  has  to  be  us«d,  as  has  been 
dene  in  this  paper. 
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Previous  work  on  the  subject  of  grid  generation  in  surfaces 
has  been  done  by  using  either  the  algebraic  techniques,  [10]  -  [12] 
or  using  the  PDE  approach,  [13]  -  [16]*.  All  these  methods  depend 
very  heavily  on  the  use  of  highly  accurate  interpolating  schemes.  On 
the  other  hand,  the  method  proposed  here  depends  only  on  the 
availability  of  the  surface  equation  in  the  cartesian  form  and  on  the 
prescription  of  the  data  on  the  bounding  curves  in  the  surface  which 
eventually  form  the  Dirichlet  boundary  conditions  for  the  proposed 
equations. 

Numerical  solutions  of  the  proposed  equations  for  the 
coordinates  in  either  simply  or  doubly-connected  regions  of  some 
known  surfaces  have  been  obtained  and  shown  in  Figs.  (1)  -  (12).  It 
has  also  been  shown  that  any  desired  control  on  the  distribution  of 
grid  spacing  can  be  excerised  by  a  proper  choice  of  the  control 
functions,  cf.  Figs.  (7)  -  (9).  Extension  cf  the  proposed  method  to 
arbitrary  surfaces  is  purely  formal. 


*  It  has  been 
be  directly 
IV. 


shown  in  [16]  that  the  equations  proposed  in  [13]  can 
obtained  by  using  equations  (4.10)  -  (4.12)  cf  sect. 
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II.  NOMENCLATURE 


b  0  »  n^.r  coefficients  of  the  second  fundamental  form 

up  ~  *tO(S 

in  the  surface  v  -  const. 

D  second-order  differential  operator;  (Eq.  3*3) 


g  * 

G  - 

v 


aet  («ij> 

6ao866  *  (go6^  1 


V 


2.  3. 


g,.  or  g  .  covariant  metric  coefficients 
ij  olp 

i  i  ri  ft 

g  J  or  g  Contravariant  metric  coefficients 

J  *  Jg  Surface  Jacobian 
v  v 

k(v),  Principal  curvatures  at  a  point  in  the 

surface  v  »  const. 

L  Second-order  differential  operator;  (Eq.  4.13) 

Unit  normal  vector  on  the  surface  v  -  const. 

P,  Q.  P“  Control  functions 

P  T 

x*  3D  Curvilinear  Coordinates 
x°  2D  Curvilinear  Coordinates 

3D  Rectangular  Cartesian  coordinates  x1  -  x,  x2  *  y.  x^  *  x 

Rectangular  components  of  n^;  xjv^  -  X^v\ 

Y(v)  v(v)  Y(v)  _(v) 

X_  «  Y  ,  X_  *  Z 
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In  this  paper  the  summation  convention  on  repeated  indices  is 

implied  when  the  same  index  appears  both  as  a  lower  and  as  an 

upper  index.  Thus  the  summation  convention  is  implied  in  T® 

but  not  in  T  .  The  summation  convention  is  also  suspended 
-  act 

when  one  repeated  index  is  enclosed  in  the  parentheses,  e.g., 

as  in  T(ct)* 
a 
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III.  THE  MATHEMATICAL  MODEL 


The  mathematical  development  of  the  model  equations  used  in 
this  paper  has  already  been  published  in  [5],  [6],  [73.  However,  for 
the  sake  of  clarity  of  exposition  we  list  here  only  the  core  steps 
leading  to  the  final  form  of  the  equations.  For  the  ensuing 
development  we  shall  constantly  use  the  conventions  and  symbols  as 
stated  in  Section  II  of  the  paper. 

The  formulae  of  Gauss,  (cf.  [4],  [17])  for  a  surface  v  - 


const,  are  written  compactly  as 


r  »  r  +  v^b 
af3  1a6-,6  -  aS 


(3.1) 


Inner  multiplication  of  (3.0  by  G 


then  yields 


Dr  +G„U<v) 


(v)_ 

[  .=  n  H, 

» 6 


(3.2) 


(3.3) 

(3.4) 


Equation  (3*2)  provides  three  second-order  partial 
differential  equations  for  the  determination  of  the  Cartesian 
coordinates  .x^,  x^  or  x,  y,  z.  However,  before  we  impose  any 
restriction  on  the  Beltramians  A^x5  for  the  purpose  of  coordinate 
control,  it  is  instructive  to  state  the  following  results: 
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(3.7) 


which  is  the  intrinsic  form  of  the  mean  curvature.  To  have  an 


extrinsic  form,  we  consider  the  formula  of  the  second  derivative  of 
the  position  vector  r  in  a  3D  Euclidean  space  given  as,  e.g.,  C7], 
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(3.8) 


=.U  '  rUc.k 


For  the  surface  xv  -  const.,  Eq.  (3.8)  is  rewritten  as 


c,<*  *  ra*  C,k' 


(3.9) 


v  k 

where  the  derivatives  with  respect  to  x  appearing  in  r  .  are  assumed 

ao 

to  have  been  evaluated  at  xv  -  const.  Taking  the  dot  product  of 

(3-9)  with  n^  and  observing  that  n^  is  orthogonal  to  any  two 

vectors  among  r  ,we  get 

n.  ( v)  rv  ,(v) 

5  •  -,afi  *  ra6  *  (3.10) 


where 


.(v)  (v)  „ 

A  «  n  .  r 

~.v 


Thus,  the  extrinsic  form  of  (3-7)  is 


k(v)  ♦  k(v)  *  eoBrv  A(v) 
I  KII  8  1 aB  A 


(3.21) 


By  using  the  form  (3.11)  in  Eq.  (3-2),  we  get  the  model 
equations  for  the  generation  of  another  surface  from  the  data  of  x, 
y.  z  ,  x,^,  y • v>  z , v  of  a  given  surface.  This  scheme  eventually  forms 
a  method  for  3D  coordinate  generation  between  any  two  (or  more)  given 
surfaces.  For  more  details  and  computational  results  on  this  aspect 
of  the  use  of  Eq.  (3.2)  refer  to  [7],  [8],  [9],  [16]. 

Since  the  subject  matter  of  the  present  paper  is  confined  to 
the  problem  of  generation  of  surface  coordinates  in  a  given  surface, 
it  is  sufficient  for  us  to  keep  the  intrinsic  form  (3-7)  in  Eq. 


(3-2) 
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IV.  THE  GENERATING  SYSTEM 


The  model  equation  for  generating  the  coordinates  in  a  given 

surface  is  now  taken  as  the  equation  (3-2).  To  be  specific,  we  take 
3 

the  surface  x  -  const,  (i.e.,  v  -  3)  as  the  given  surface.  In 

12  3 

expanded  form,  taking  r  -  (x.y.z),  x  -  x  -  n,  x  -  ?,  the  three 
equations  are 


Dx  ♦  t(43)£)  x  ♦  (43)n)  x  ]G.  -  X^3)R, 
*  t  «  n  3 

Dy  *  C(a(3)e)  y?  +  (43Ji  )yn3G3  -  y(3)r, 


Dz  ♦  C(a^3)e)  z_  ♦  (A*3)n)  z  ]  G,  -  z(3)r, 
<=  t  *  n  3 


(4.3) 


where 


.(3)  r  ^  /Jr  1  ^  //g, 

l2  a  ar  12  3  3 

~^2g12^12  ”  g1ir22  “  g22Tll)/G3’ 

(4.6) 

We  now  suitably  restrict  the  distribution  of  the  Beltramians 
so  as  to  have  available  a  sort  of  flexibility  in  the  choice  of  the 
coordinates  £,  n  in  the  surface.  The  most  general  form  which  can  be 
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-  1  - V  \  1  V*  V 


chosen  (cf.  [6],  is 


g 

g 


agp1 

a8* 


a8  p2 
a8’ 


(4.7a) 

(4.7b) 


1  2 

where  PQg,  p^g  are  symmetric  in  a,  8  and  represent  six  arbitrary 
chosen  control  functions.  Using  the  summation  convention  in  (4.7) 
and  the  formulae 


g11  -  822/G3<  g12 


2^ ^ 9  ^ 


22 


g11  /C3 


we  have 


a(3). 
a2  5 


r  A(3)n 

V  2  n 


(4.8) 


where 

p  “  g22P 11  ~  2gl2P12  +  g11P22’  (4.9a) 

Q  *  g22P1 1  ~  2g12P12  +  gllP22‘  (4.9b) 


The  final  form  of  the  eqs.  (4.1)  -  (4.3)  is  now 


Lx  -  X(3)R  , 

(4.10) 

Ly  -  Y^3)R  , 

(4.11) 

Lz  =  Z(3)R  , 

(4.12) 

where 

1  "  g223«  '  2g123^n  +  9lldn»?  + 


G3C(A23>t)36  +  <43>n)V 


(4.13a; 


-  e22s„  -  2g12as„  *  g„»nn  *  p»{  *  (4  Ub 

Equations  (4.10)  -  (4.13)  are  the  basic  generating  equations 
for  the  curvilinear  coordinates  in  a  given  surface.  The  function  R 
appearing  on  the  right  hand  side  of  the  Eqs.  (4.10)-  (4.12)  is  given 
by 

R  -  (kj3)  +  kj3))G3- 

The  quantity  in  the  parentheses  is  twice  the  mean  curvature  of  the 
surface  and  as  noted  in  Sect.  Ill  is  invariant  to  the  coordinates 
introduced  is  the  given  surface.  Consequently  k.£3^  ♦ 
reflects  a  basic  geometrical  aspect  of  the  surface  and  is  a  function 
of  the  coordinates  x,  y,  z.  Since  the  Cartesian  form  of  the  equation 
z  *  f  (x,y)  is  assumed  to  be  available  for  the  surface  under 
consideration,  it  is  obvious  from  elementary  differential  geometry 
that 


x(v)  ♦  k(v) 
I  II 


[  ( 1  ♦  z  )  z  -2z  zz  + 
y  xx  x  y  xy 


(1  ♦  z2)  z  ]/(1  +  z2  +  Z2)3/2 
x  yy  x  y 


(4.14) 


For  arbitrary  surfaces,  it  is  always  possible  to  generate  the 
Cartesian  form  z  *  f(x, y)  of  the  surface  by  least  square  data  fitting 
thus  having  k^  ♦  k  as  a  function  of  x  and  y.  It  is  also  possible 
to  solve  only  Eqs.  (4.10)  and  (4.11)  while  calculating  z  from  the 
given  equation  z  -  f(x,y)  of  the  surface.  All  these  aspects  have 
been  discussed  in  Sect.  VI. 
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V.  GENERATING  EQUATIONS  BASED  ON  THE  VARIATIONAL  PRINCIPLE 
Two  attractive  features  of  the  chosen  equations  (4.10) 
(4.13)  are  their  simplicity  and  their  explicit  dependence  on  a  basic 
surface-  geometric  property,  viz.,  the  mean  curvature.  Any  other  ad 
hoc,  though  consistent,  set  of  equations  can  also  be  used  to  generate 
the  coordinates  but  then  it  will  be  extremely  difficult  to  isolate 
those  basic  surface-geometric  properties  which  distinguish  one 
surface  from  the  other.  In  this  connection  the  variational  approach 
is  also  a  possibility  which  has  been  used  by  Brackbill  et  al  [19]  for 
coordinate  generation  in  a  flat  space,  i.e.,  a  2D  plane. 

Generally,  let  us  consider  the  surface  functional 

1  (5.1) 

Then  51  =  0  leads  one  to  the  Euler-Lagrar.ge  equations 


a*r  i 


3(/G^<£  ^ 


(5.2) 


where,  referring  to  Sect.  II  we  have  used  the  summation 

I 

x  (r  =  1,2,3)  as  the  rectanguler  Cartesian  coordinates, 
r 


r,B 


a* 

3x 


_r 

0’ 


3% 


r,a6 


3xa3x6 


convention, 

and 


With  these  notations,  it  is  a  direct  algebraic  problem  to  show  that 

/G  8  i  9G3 

Lx  =  — —  — —  (—  , - ±  )  , 

r  2  3xT  /g  3xr  v 

3X  °3  r,Y  (5.3) 


where  the  operator  L  has  earlier  been  defined  in  (4.13a).  Let  i  be  a 
function  only  of  x  ,  then  expanding  (5.2)  while  using  (5.3),  we  get 
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(5.4) 


where 


Lx 


-M 


M  -  ± 
2 


.Jh  _U 

3xr,s  ax6 


a4> 


3G. 


3x 


r,B  3x 


)+Z?±J*±  \ 

<t>  d*.  .  ' 

>,p 


8  <P  dx 


(5.5) 


Thus  the  generating  system  (5.4)  Is  of  the  same  form  as  originally 
proposed  but  it  looks  to  be  a  formidable  problem  to  select  that  4> 
which  yields  the  right  hand  sides  of  Eqs.  (4 . 10)-(4. 1 3) .  One  simple 
result  for  the  case  when  <p  *  1  is  obvious.  For,  in  this  case  the 
minimization  of  I  in  (5.1)  implies 


and  these  are  the  equations  for  a  minimal  surface.  Another  case  in 
which  <p  =  F/Cj  with  F  still  as  a  function  only  of  x  yields  the 
Euier-Lagrange  equations  as  (5-4)  with  M  defined  as 


M 


1  .  3G3 

2F  1 3x 


3F 


3G. 


3F , 


r,B  3x 


B 


3x8  3xr,B 


dF 


1  3G3  fi_3  f _ 

2G3  3xr, B  3x8  F  3xS  3x 


r,B 


(5.6) 


The  above  case  of  <p  -  F/G^  is  of  interest  because  the  choice 

F  “  g1 1  +  g22’ 
or 


11  22 
4>  -  g  ♦  g 


(5.7a) 

(5.7b) 


is  equivalent  to  the  "smoothness"  problem  in  2D  plane  coordinates  as 


shown  by  Brackbill  [19]. 

It  must, 

however, 

be 

stated 

that 

"smoothness"  of  coordinates  in 

a  2D  plane 

problem 

is 

due  to 

the 
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VI.  NUMERICAL  IMPLEMENTATION 


The  numerical  solution  of  Eqs.  (4. 10)-(4. 12)  can  be  obtained 
by  any  suitable  numerical  method  of  solution  which  has  proved  useful 
in  any  elliptic  grid  generation  problem.  In  this  paper  the  equations 
have  been  discritized  by  using  central  differences  for  both  the  first 
and  second  derivatives  and  then  solved  iteratively  from  an  initial 
guess  by  using  the  LSOR.  The  main  difference  between  the  coordinates 
in  a  flat  space  and  in  a  surface  is  the  appearance  of  the  right  hand 
side  terms  in  which  the  quantity  R  can  be  established  a'  rpriori. 
This  requires  a  knowledge  of  the  equation  of  the  surface  z  -  f(x,y), 
which  when  used  in  (4.14)  yields  R  as  a  function  of  x  and  y.  For 
arbitrary  surfaces  the  equation  z  =  f(x,y)  can  be  established  by  the 
least  square  method,  [18]. 

To  demonstrate  the  potential  of  Eqs.  (4 . l0)-(4. 12)  as  a  viable 
set  of  equations  for  the  generation  of  surface  coordinates,  we  have 
selected  three  well  known  surfaces  for  the  purpose  of  introducing  a 
desired  system  of  coordinates  in  them. 

a.  Coordinates  In  An  Elliptic  Cylinder  Forming  A 
Simply-connected  Domain. 

This  problem  is  a  prototype  of  coordinate  generation  in  a 

given  piece  of  a  surface.  The  region  under  consideration  forms  a 

simply-connected  region  bounded  by  the  space  arc  n  =  nQ,n  =  n1 ,  £ 

£  ,  and  £  *  5, .  Here  n  =n  , n.  are  the  elliptical  arcs  in  the 

o  1  0  1 

xy-plane,  and  £  =  £  are  straight-lines  parallel  to  the  z-axis. 
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and  5  =  •  Here  n  =  no»n.|  are  the  elliptical  arcs  in  the 

xy-plane,  and  £  =  £  are  straight-lines  parallel  to  the  z-axis. 


The  equations  are: 


C  =  5o  :  x  =  -a.  y  =  0,  Zj  4  z  4  ze , 

5=6,  :  x  »  a,  y  =  C,  Z]  <  z  <  zo  . 

The  Dirichlet  boundary  conditions  are  provided  by  the  data  of  (6.1) 
for  the  solution  of  Eqs.  ( 4 . 1 C) -( 4 . 1 2 ) . 

b.  Coordinates  In  An  Arbitrary  Ellipsoid  Cut  By  Tne  Planes  z  =  z  , 
z  =  z,. 

This  case  is  of  coordinate  generation  in  a  dcubly-connrcted 
region  bounded  by  two  closed  space  curves  on  an  ellipsoid.  Tne  space 
curves  r,  *  nQ  and  n  *  ,  are  giver,  by 


2  ? 

n-n.  :x_*z2+z__=l,z  =  z.. 

a2  b2  c2  (6.2) 

We  now  imagine  a  cut  joining  the  curves  n  *r,  ,  n  =  n1  whil*-  still 

remaining  in  the  surface.  As  in  the  CD  case  no  boundary  cord :  -  icrs 

can  be  prescribed  on  the  cut  line.  However,  sine*  the  v  :)u-  s  of 

x,y,z  above  and  below  the  cut  should  be  th'  sa~o,  impose  tne 

periodicity  conditions: 
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x^.n)  -  x(£o,n).  y(Crn)  -  yU0,n).  zU1f  n)  -  zUo,n).  (6.3) 

The  Dirlchlet  conditions  (6.2)  and  the  periodicity  conditions  (6.3) 
yield  a  unique  solution  to  the  set  of  equations  (4.10)  -  (4.12). 

C.  Coordinates  In  An  Arbitrary  Elliptic  Parabolid  Cut  By  The  Plane 

2  “  V  z  -  V 

This  is  again  the  case  of  coordinate  generation  in  a 
doubly-connected  region  bounded  by  two  closed  space  curves  on  an 
elliptic  parabolid.  The  space  curves  n  *  nQ  and  n  *  n1  ,  are  given  by 

n  -  n  :  ~  ♦  4  =  z  ’  Z  "  2°’ 

o  a  .2 
b 


n  -  n. 


2  2 
2L_  ♦  Z_  =z 
2  .2 
a  b 


(6.4) 


Under  the  boundary  data  (6.4)  along  with  the  periodicity  conditions 

(6.3)  the  equations  (4.10)  -  (4.12)  have  been  solved. 

In  all  cases  (a)  -  (c)  the  control  functions  P1 „  and  p2 

aB  aB 

(3)  (3) 

have  been  set  equal  to  zero,  i.e.,  A  2  n  -  0,  &2  n  *  0,  and  the 
results  are  demonstrated  in  Figs.  1  -  6  and  Figs.  10  -  12.  Figures 
7-9  show  the  results  of  a  coordinate  concentration  near  the  curve  z 
*  zq  *  0.9  of  case  (b).  In  this  case  we  have  taken 


2  2 
°*  P11  "  P?2 
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where  *  -  1.1  Is  a  constant. 

The  computer  programs  which  have  been  developed  to  solve  the 
equations  (4.10)  -  (4.12)  have  been  used  successfully  for  all  the 
cases  enumerated  above  both  with  and  without  coordinate  contraction. 
Also  each  case  was  repeatd  to  determine  whether  it  is  necessary  to 
solve  the  z-equation  too  along  with  the  x  and  y-equation.  It  has 
been  found  that  solving  all  the  three  equations  (4.10)  -  (4.12)  or 
solving  only  (4.10)  and  (4.11)  while  taking  the  z  iteratively  from  z 
*  f(x,y)  produces  almost  the  same  results. 


V.  Conclusions 


A  set  of  second  order  partial  differential  equations  have  been 
developed  and  then  solved  numericaly  to  generate  the  coordinates  in  a 
given  surface.  The  proposed  equations  are  a  logical  outcome  of  the 
formulae  of  Gauss  and  thus  explicity  depend  on  some  basic 
differential-geometric  properties  of  the  surface  in  which  the 
coordinates  are  to  be  introduced.  The  proposed  method  of  surface 
coordinate  generation  is  simple  to  implement  and  is  capable  of 
extension  to  arbitrary  surfaces. 
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Figure  1.  A  simply-connected  region  on  the  surface  of  an  elliptic  cylinder 
data;  a  =  1.0,  b  =  0.5,  and  z  =  zq  =  2.0,  z  =  •  0.0.  Viewed  after  a  20* 
clockwise  rotation  about  the  x-axis. 


1  jure  3 


Data  same  as  in  Fig.  1. 


Viewed  after  a  20"  clockwise  rota 


ion  about 


e  :-axis. 


Figure  4.  A  doubly-connected  region  on  the  surface  of  an  ellipsoid 
a  *  5.0,  b  =  3.0,  c  =  1.0,  and  cut  by  the  planes  z  =  zQ  -  0.9,  z  = 
Viewed  after  a  20 0  clockwise  rotation  abcu  the  <-axis. 


axis 


4,  coordinate  contraction 


Figure  10.  A  doubly-connected  region  on  the  surface  of  an  elliptic  paraboloid 
data:  a  =  2.0,  b  =  1.0,  and  cut  by  the  planes  z  =  zQ  =  0.01,  z  =  z^  *  1.96. 
Viewed  after  a  20;  rotation  about  the  x-axis. 
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ABSTRACT 


This  paper  Is  devoted  to  the  numerical  solution 
of  a  set  of  second  order  elliptic  partial  differential 
equations  for  the  generation  of  three-dimensional 
curvilinear  coordinates  between  arbitrary  shaped  bod¬ 
ies.  Starting  from  the  given  data  on  the  inner  and 
and  outer  surfaces,  a  series  of  surfaces  are  generated 
which  are  connected  with  one  another  in  such  a  manner 
chat  a  sufficiently  smooth  three-dimensional  coordi¬ 
nate  net  is  realized.  A  cumber  of  cases  pertaining 
to  the  numerical  grid  generation  between  two  given 
surfaces  have  been  discussed. 

INTRODUCTION 

Grid  generation  by  computer  simulation  through 
carefully  selected  mathematical  models  has  become  a 
useful  tool  in  many  branches  of  engineering  and  phy¬ 
sics.  Although  the  developments  In  grid  generation 
are  mostly  due  to  the  problems  In  computational  fluid 
dynamics,  the  techniques  developed  sre  also  appli¬ 
cable  to  all  those  areas  where  field  equations  sre 
tu  be  solved  on  complicated  domains.  The  main  idea 
behind  grid  generation  by  computer  simulation  is  to 
have  a  well  structured  mathematical  model  which  is 
solved  for  the  physical  coordinates  (the  Cartesian 
coordinates)  as  functions  of  the  curvilinear  coordi¬ 
nates  under  the  constraints  of  the  specified  bound¬ 
aries  on  which  the  physical  coordinates  are  known  in 
advance. 

In  tiie  generation  of  two-  and  three-dimensional 
coordinates  around  and  between  bodies  of  arbitrary 
shape  there  exists  two  main  approaches:  (1)  alge¬ 
braic  methods  methods  based  on  solutions  of 
partial  differential  equations  (Preferably  elliptic 
partial  differential  equations).  The  present  method 
is  concerned  with  the  second  class,  that  is,  on  solu¬ 
tions  of  elliptic  partial  differential  equations. 

Refer  to  Thompson,  Wars!,  and  Mas tin  (1)  for  more 
information  on  algebraic  methods. 

Elliptic  grid  generation  systems  have  some  ad¬ 
vantages  which  seem  to  make  them  a  better  approach  to 
take  (2).  The  resulting  grid  is  naturally  smoother 
and  there  is  less  danger  of  overlapping  of  the  grid 
lines.  Because  of  the  elliptic  nature  of  the  equa¬ 
tions  boundary  slope  discontinuities  do  not  propagate 
into  the  field.  Grids  generated  by  solutions  of 
elliptic  equations  are  relatively  easy  to  adapt  to 
gtr., ral  boundary  configurations  and  they  have  the 
ability  to  incorporate  features  such  as  concentration 
-f  arid  lines,  smoothness,  and  orthogonality. 

Tbit  work  is  devoted  to  the  generation  of  three- 
dittrsicnal  grids  between  two  arbitrary  shaped 
bodies.  The  method  has  its  origin  in  the  formulae  of 
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Gauss  for  a  surface.  In  this  scheme  an  inner  and 
outer  boundary  surface  Is  assumed  to  be  Initially 
given.  The  inner  bodies  investigated  are  those  of 
various  sites  of  ellipsoids,  fuselage  shapes,  two 
intersecting  spheres  or  ellipsoids,  and  a  wing-body 
combination,  la  all  cases  the  outer  surface  is 
either  a  sphere  or  an  ellipsoid.  A  series  of  sur¬ 
faces  are  generated  between  the  inner  and  outer 
boundary  surfaces  in  such  a  way  that  only  two  coordi¬ 
nates  vary  in  each  surface.  The  equations  solved 
(Eqs.  (4))  are  structured  in  such  a  way  that  there  is 
an  automatic  connection  between  succeeding  generated 
surfaces.  The  result  is  a  sufficiently  differenti¬ 
able  three-dimensional  coordinate  system  in  the  en¬ 
closed  region.  Provisions  are  made  for  contraction 
of  coordinates  in  the  generated  surfaces  and  results 
are  shown  for  contraction  near  the  inner  body. 

In  the  case  of  a  wing-body  combination  being  the 
inner  body  a  set  of  partial  differential  equations 
(Eq.  (6))  are  solved  which  transform  whatever  coordinate 
system  la  initially  given  on  the  body  to  one  In  which 
one  of  the  coordinate  lines  Is  the  trace  of  the  wing 
and  body  intersection.  This  new  coordinate  system 
allows  for  a  smooth  matching  of  the  wing  with  the 
body.  With  this  now  as  the  coordinate  system  on  the 
inner  body  the  equations  are  solved  to  generate  the 
three-dimensional  coordinate  system  between  the  wing- 
body  and  the  outer  boundary  surface. 

MATHEMATICAL  MODEL 

The  basic  mathematical  development  of  the  method 
to  be  discussed  is  available  in  Warsl  (3),  (4),  and 
the  details  on  the  computer  simulation  are  available 
in  Ziebarth  (5),  and  Tlam  (6).  Below  we  shall 
briefly  state  the  essential  structure  of  the  proposed 
model. 

As  mentioned  in  the  introduction,  the  present 
model  has  its  origin  In  the  formulae  of  Gauss  for  a 
surface.  The  formulae  of  Gauss  are  essentially  a  sat 
of  compatibility  relations  connecting  the  second 
derivatives  with  the  surface  Christoffel  symbols,  the 
coefficients  of  the  second  fundaaentsl  form,  and  the 
unit  nOTtal ,  all  evaluated  in  the  surface.  Denoting 
the  three  Cartesian  coordinates  in  space  by  r  ■  (x, 
y,z)  and  the  curvilinear  coordinates  by  x*  * 

(i  •  1,2,3),  the  formulae  of  Gauss  for  the  surface 
xv  •  const,  are 

!ob  ‘  Ta8  ?6  +  boe  !(v)  •  <l) 

where  repeated  lower  and  upper  indices  Imply  summa¬ 
tion,  and  , 


Rx0*8  *  -*  ’  Rx4 
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In  Eq.  (X)  T*.  are  che  surface  Chriscoffel  symbols 
defined  aa  08 


ix6 


equations  for  the  generation  of  surface  coi-rdinates 

are 

D!  +  P*22fi  +  QSjlfn  ’  ?R  •  (4> 


where  g„g  are  the  coefficients  of  the  first  fundamental 
f orn .  The  coefficients  of  che  second  fundamental 
form,  bQg  ,  are  defined  as 


( v)  V 

where  the  unit  normal  n  on  che  surface  x  •  const. 


(!o  *  x  fs1  • 


Referring  to  Fig.  1,  let  the  given  inner  and  outer 
body  surfaces  be  given  as  n  *  hg  and  ^  on  which 
the  values  of  r  *  (x,y,r>  are  known.  Since  the  coor¬ 
dinates  •;  and  ;  have  been  chosen  in  these  surfaces, 
the  derivatives  r.  on  the  inner  and  outer  surfaces, 
denoted  as  (r  )g‘and  (r^),,  are  known.  To  connect 
one  generated^surface  through  Eq.  (4)  with  the  other 
we  specify 

U  ’  fl(r)(ft)B  +  f2(r)(f;>-  <5) 


It  oust  be  mentioned  chat  che  range  of  the  Creek 
indices  ia  only  on  two  values.  Thus,  for  v  ■  3,  the 
index  range  is  from  1  to  2.  From  (1),  writing  the 
equations  for  a  •  1,  8*1,  o  -  1,  8*2,  and  a  ■  2, 
6*2,  and  multiplying  Che  equacions  respectively  by 
*22’  *12  *nd  >21  *nd  a<ldin8  them  together,  we  get 


where  fjf")  and  f^Cr.)  are  suitable  weights  which  are 
such  that 


Dr  ♦  [d2f)r{  +  (Ajitfr  JCj  -  nR,  (2) 

where  chese  equations  are  for  the  surface  (  •  const, 
on  which  £  and  n  are  the  current  coordinates.  In 
Eq.  (2)  a  variable  subscript  Implies  a  partial  deriv¬ 
ative,  and 

D  “  *223ti  ‘  2*i23£n  +  *ll3nn  * 

’  T^p  lH  (*223c  ‘  *123n^ 

/C3  /C3 

+  **  J~cT  *  *12  V1  ’ 

V  C3 

C3  “  *11*22  -  (*12)2  • 


n  -  (X.T.Z)  , 

Fig.  1  (a)  Topology  of  the  given  surfaces, 

(b)  Surface  to  be  generated. 


*  *  (kj  +  *pC3  • 

where  kj  ♦  lu  is  twice  the  mean  curvature  of  the 
surface  at  any  point.  The  quantities  and  Ago 
ere  the  surface  Laplacians  or  Beltraaians. 

A  sec  of  generating  equations  can  now  be  obtained 
by  setting 

V  “  p«22/c3  .  (3) 

42n  "  Q*ll/C3  ’ 

where  P  and  Q  are  arbitrary  specified  functions  which 
exert  the  desired  control  on  the  distribution  of 
coordinates  in  the  surface  to  be  generated.  Substi¬ 
tuting  (3)  in  (2),  the  elliptic  partial  differential 


The  set  of  Eqs.  (4)  essentially  complete  the 
development  of  a  consistent  mathematical  model  for 
the  generation  of  surface  coordinates,  where  the 
means  of  proceedings  from  one  generated  surface  to 
the  other  is  provided  by  Eqs.  (5). 

SURFACE  COORDINATES 

Equations  (4)  can  also  be  used  in  a  straight 
forward  manner  to  yield  surface  coordinates  based  on 
prescribed  data  on  certain  contours  of  the  surface 
Itself.  These  equations,  (7),  provide  a  means  of 
introducing  a  new  coordinate  system  based  on  the 
data  of  an  existing  coordinate  system,  e.g.,  from  a 
polar  coordinate  system  of  Craidon  (8)  to  an  0-type 
coordinate  system  in  the  surface,  which  may  be  needed 
In  certain  wing-body  combinations.  These  equations 
are  stated  as  follows. 


SB 


“«  -  2buU  +  cu«  ♦  '“t  +  *  J2V  • 
*vu  '  2bvu  +  CVCi  +  *1  +  ‘  J2V  * 


(6.) 


where  u,  v  are  the  parametric  coordinate*  In  the  sur¬ 
face  which  are  to  be  transformed  to  the  new  coordi¬ 
nates  (C.O,  and  the  coefficients  a,  b,  etc.,  are 

•  ’  <iuv^  ♦  2g13«tvt  ♦  I33u{)/52  . 

b  •  <»llv;vt  ♦  *13(\v(  +  “tV  +  *33Vc,/52’ 

c  -  (iuv?  ♦  2il3Vc  +  i33u2)/C2  * 


•  G, 


t  C 


‘  Vc 


(6b) 


and  i2u,  i,v  are  the  Beltramians  defined  earlier. 

An  overhead  bar  meanr  a  metric  coefficient  in  the  u, 
v  system,  e.g. 

2.22 

g, ,  •  x  ♦  y  +  *  etc. 

"11  u  'u  u 


Equations (6a)  and  (6b)  have  been  programmed  in 
Refs.  (5)  and  (6).  In  (S)  these  equations  have  been 
used  to  provide  those  coordinates  which  are  such 
that  the  inner  0-type  coordinate  conforms  to  the  con¬ 
tour  created  by  the  intersection  of  a  wing  with  the 
fuselage.  Equations  similar  to  Eqs-  (6a)  and  (6b)  have 
been  derived  by  Caron  and  Camarero  (9)  and  Whitney 
and  Ihomas  (10),  the  later  reference  has  followed  a 
very  teoius  procedure. 

NUMERICAL  SOLUTION  OF  3-0  EQUATIONS 

The  system  of  partial  differential  equations, 

(Eqs.  (4),  for  the  generation  of  a  three-dimensional 
coordinate  system  between  two  arbitrary  shaped  bodies 
has  been  solved  by  finite  differences  using  point- 
SOR.  Initially  the  coordinates  on  the  two  boundary 
surfaces  must  be  prescribed  by  some  method.  In  this 
research  these  coordinates  were  prescribed  analyti¬ 
cally  for  simple  shapes  or  were  calculated  by  a  com¬ 
puter  subroutine  of  Craldon  (8).  Although  the  choice 
of  an  outer  boundary  surface  is  arbitrary  there  must 
be  a  one-to-one  correspondence  between  points  on  the 
two  boundary  surfaces.  In  addition,  these  points 
need  to  be  chosen  so  that  the  lines  Joining  the 
points  cn  the  inner  and  outer  surfaces  do  not  Inter¬ 
sect. 

Regardless  of  the  method  chosen  to  correspond 
points  between  the  inner  and  outer  boundary  surfaces 
(see  (5)  for  a  discussion  of  these)  it  has  been 
found  to  be  a  shape  dependent  mechanical  exercise. 

For  example,  Figures  2  and  3  show  the  correspondence 
between  points  on  the  Inner  and  outer  surface  for 
the  case  of  Intersecting  ellipsoids.  Although  in 
both  cases  there  is  a  one-to-one  correspondence  be¬ 
tween  points  cn  each  surface,  when  the  correspondence 
of  Fig.  2  is  used  the  solution  will  not  converge, 
while  that  of  Fig.  3  does  yield  a  solution  to  the 
differential  equations.  There  are  two  main  problems 
with  the  correspondence  of  Fig.  2.  The  lines  near 
the  right  end  are  close  together  and  appear  to  be 
parallel.  This  is  believed  to  be  the  reason  for  the 
divergence  uf  the  solution  since  the  computations 
stopped  in  this  region.  Another  problem  with  the 


Fig.  3  Correspondence  of  points  which  allows  conver¬ 
gence  of  the  solution. 

correspondence  in  Fig.  2  is  that  the  correspondence 
line  labeled  A  in  Fig.  2  is  not  vertical.  This  line 
should  be  chosen  to  be  along  the  y-axis  (x  -  z  •  0) 
so  that  it  is  not  permitted  to  change  from  a  quadrant 
where  x  is  negative  to  one  where  x  is  positive.  This 
second  correspondence  problem  was  found  to  cause  even 
the  solution,  for  a  simple  ellipsoid  as  the  inner 
body,  to  diverge. 

Runs  were  made  with  both  an  ellipsoid  and  a 
sphere  as  the  outer  body.  Ko  advantage  was  found  In 
taking  the  outer  body  to  be  an  ellipsoid  so  a  sphere 
was  used  In  all  other  cases. 

In  conclusion,  the  correspondence  between  the 
inner  and  outer  body  is  important  for  the  solution  of 
equations  (6)  to  converge.  There  does  not  seem  to  be 
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a  general  method  which  works  (or  all  inner  body 
shapes,  therefore,  the  zest  useful  way  seems  to  be 
a  geometric  manipulation  which  is  dependent  on  the 
inner  body  configuration. 


Once  the  correspondence  has  been  established 
between  the  inner  and  outer  bodies  an  initial  guess 
must  be  made  for  the  values  of  the  Cartesian  coordi¬ 
nates  in  the  field.  In  all  cases  linear  Interpola¬ 
tion  between  the  points  on  the  inner  and  outer  sur¬ 
faces  was  used  successfully  to  accomplish  this. 


where  <  *  1.0  implies  no  contraction  : -or  the  inner 
body  and  •:  slightly  greater  than  1.9  (*.g.  1.05  or 
1.10)  implies  contraction.  Here  ”  is  e.nu-e.-a :cj  as 
ar.  integer. 

In  the  case  of  mere,  complicated  shapes  such  as 
intersecting  ellipsoids  the  exact  solution  coy  not  b« 
kr.mwn  between  it  and  an  outer  surface.  However, 
fer  certain  lines  in  the  field  equations  (-)  can  be 
solved  in  their  limiting  ferms.  One  s-ch  lire  is 
the  line  A'B'  in  Fig.  5.  Wars!  (■*)  has  show-  thet 
along  this  line  the  expression  for  x  can  be  written 


t  or  K  /  ‘  “A  t  -  0  or  X  .  1 


Fig.  4  Lines  c  «  1  and  ;  -  ;max  for  a  prolate 
ellipsoid  surrounded  by  a  sphere. 

Before  a  solution  algorithm  can  be  applied  to 
equations  (4)  appropriate  boundary  conditions  for 
i  »  1  and  ;  *  (Fig.  4)  must  be  determined.  For 

the  case  of  a  prolate  ellipsoid  surrounded  by  a 
sphere  an  exact  solution  exists  and  is  used  to  cal¬ 
culate  these  boundary  conditions.  The  solution  is 
given  by,  (4), 

*  "  [A  e®n^  ♦  C]cos;  (7) 


Fig.  5  Intersecting  ellipsoids  s-r receded  by  a 

sphere. 


x(n)  ■  xB  + 


At51’’ 


(e  -  cosh-,a)sinhrig 


e  -  5inhn0 


r)m  _ 1 _ 

B  -  tn{J— - !r— ■]*-  '  \  ) 
‘sinhn.  B 


e  (coshhg  -  sinhbg) 


•  -  slnhn. 


D  ■  sinhnB 

and  n(n)  may  be  taken  as 


and  the  subscripts  B  and  *>  refer  to  the  inner  and  j 
outer  body  respectively.  The  expression  for  x  in  | 
(8)  has  been  used  in  the  case  of  intersecting  sphere^ 
and  ellipsoids  to  evaluate  the  boundary  conditions  : 
on  the  right  segment  (AA'ln  Figure  5).  A  value  of  <3 
slightly  larger  than  1.0  is  required  in  (V)  to  pro-  1 
duce  coordinate  systems  which  are  smooth  near  the  .) 
right  segment.  Because  of  the  orientation  of  the 
intersecting  bodies  with  the  Cartesian  coordinates 
the  boundary  conditions  for  the  left  segment  (3B’ 
in  Figure  5)  can  be  determined  by  an  expression  like* 

(7). 

Once  the  coordinates  on  the  inner  and  outer  1 
bodies  have  been  specified  the  t-derivatives  on  thesa 


60 


Of  J 


tur(’.as  can  also  be  calculated  by  second  order  finite 
differences.  These  need  to  be  calculated  on  both  the 
inner  anc  outer  surfaces  for  all  values  of  K  from  K 
1  .  ;  to  K  *  UMAX-1  and  for  all  values  of  1  from  1*1 

t0  ;  ■  1MAX  (Fig.  b) .  The  (-derivatives  are  not 
retired  at  K  •  1  and  K  •  UMAX  since  these  are  bound¬ 
er*  conditions. 

Special  care  Bust  be  taken  at  locations  where 
there  exists  a  mathemat i cal  discontinuity  (e.g.,  at 
the  intersection  point  K  *  Ke  In  Fig.  5).  A  simple 
averaging  of  the  (-derivatives  through  this  region 
wt,  found  to  adequately  smooth  out  these  derivatives. 

Two  mere  concerns  need  to  addressed  before  the 
differential  equations  can  be  solved.  The  first  is 
the  spec  if ica tion  of  the  functions  f^(r)  and  f,('i) 
appearing  in  equation  (5).  In  this  research  fj  was 
taker,  to  be  a  linear  function  of  z  where 


r,  -  n 

m 


and  f-(z)  *  1  -  fj(z).  The  second  concern  is  the 
specification  of  the  redistribution  functions  P  and 
Q  in  equation  (4).  In  this  work  only  contraction  of 
the  --luies  near  the  inner  body  was  considered, 
therefore,  cr.ly  the  function  Q  was  specified  and  P 
was  set  equal  to  zero.  Based  on  earlier  two- 
dimensional  work  (11),  (12),  Q  was  prescribed  as 

i, ,  U  +  ('  -  hj  Inline 
Q  .  - - _ -  . 

1  +  (*  -  hg) In  < 

KTbSr.lCAL  SOh'-'TION  OF  SURFACE  EQUATION’S 


Fig.  6  Nomenclature  for  ellipsoid  as  an  inner 
boundary  shape. 


Fig.  7  Coordinates  {,  (,  n  on  the  wing-body  com 
tion . 


In  the  case  of  a  wing-body  combination  equations 
(4)  are  solved,  but  only  after  a  transformation  of 
the  coordinate  system  on  the  body  is  made.  The  body 
is  initially  covered  by  a  coordinate  system  (u.v). 

By  initially  solving  equations  (6)  the  (u.v)  coordi¬ 
nate  system  is  transformed  into  a  (i,()  coordinate 
system,  on  the  same  body.  The  ((,;)  coordinate  sys¬ 
tem  will  cover  only  half  of  the  original  body,  but 
it  has  the  advantage  of  having  one  of  its  coordinate 
lines  as  the  trace  of  the  intersection  between  the 
body  and  the  wing.  Fig.  7 . 

In  order  to  solve  the  differential  equations  (6) 
boundary  conditions  need  to  be  chosen  on  the  surface 
of  the  body.  Although  the  method  for  picking  points 
is  arbitrary  they  need  to  satisfy  a  few  basic  re¬ 
quirements.  The  inner  boundary  is  to  be  the  trace 
of  the  intersection  line  between  the  body  and  wing 
(line  labeled  A  in  Fig-  8;.  The  outer  boundary 
(line  B  in  Fig.  8)  should  be  very  close  to  the  line 
which  cuts  the  body  into  two  similar  halfs  (line  C 
in  Tig.  8).  It  can  even  be  coincident  with  part  of 
this  curve  as  shown  in  Figure  8.  However,  the  outer 
boundary  should  not  be  taken  to  be  entirely  on  this 
line  since  It  needs  to  have  a  variation  over  the 
surface  of  the  body,  i.e.  ,  a  variation  in  all  three 
Cartesian  coordinates  x,  y  and  z. 

When  picking  the  points  on  the  inner  or  outer 
boundary  it  is  necessary  to  be  able  to  find  both  the 
Cartesian  coordinates  (x,y,z)  and  the  (u.v)  coordi¬ 
nates  for  points  within  any  patch.  This  is  accom¬ 
plished  by  considering  parameters  -  and  v  which  vary 
from  0  to  i  within  any  patch.  Here  w  varies  ir.  the 
direction  of  the  u  coordinate  and  v  varies  in  the 
direction  of  the  v  coordinate.  Obviously  since  u 


Fig.  8  Choice  of  boundary  curves  or.  inner  body. 


varies  from  0  to  r  then 


where 


and 


n 


X 


s 

NPU 


NPV  *  number  of  patches 
U  direction 


v 


6) 


where 


(1,0) 


NPV  •  number  of  patches  in  V 
direction 

Based  on  the  parameters  u  and  v  a  blending  equation 
can  be  defined  for  these  Surface  patches  similar  to 
the  one  used  by  Craidon  (8).  This  bicubic  surface 
patch  equation  has  the  form 

r(w,v)  *  F1(u)Fi(  .)r0  0  +  F1(u)F2(v)rQ  2 

+  F1(u)F2(u)r1  Q  ♦  F2(.)F2(.)r1  2 

♦  F,(v)F,(u)  r  +  F,(v)F,(w)r 

1  3  -“0,0  2  3  '“0,1 


+  F,(v)F,(u)r  +  F,(u)F.(i.)r 

l  “  -U10  a  «  ,u1(l 


Fig.  9  (a)  Two  of  the  surface  patches  with  labeling 

of  corner  points,  (b)  Isolated  surface 
patch  i,j  showing  parameters  .  and  .. 


+  F2(l-)F3(v)rv  +  F2(u)F4(u)rv 


+  F,(w)F,(v)r  +  F,(u)F.(v)r 

'1,0  “  1,1 


+  F,(u)F,(v)r  +  F,(u)F,(v)r 

3  3  0  i  *»  ^ 


+  F3U)F4('“)!uv1  fl  +  x 


where  the  blending  functions  can  be  taken  as 


Fx(9)  -  263  -  392  +  1 
F2(S)  -  392  -  293 
F3(8)  -  63  -  292  +  8 
F.(8)  -  83  -  82 

H  w 


and  6  Is  a  dummy  variable.  The  mixed  second  deriva¬ 
tives  of  r  were  dropped  from  equation  (12)  in  all 
calculations.  As  was  discussed  earlier  in  this 
chapter  r  and  also  the  derivatives  of  r  with  respect 
to  u  or  v  can  be  calculated  at  each  corner  by  finite 
differences.  Thus,  equations  (14)  can  be  used  to 
find  f  at  any  point  in  any  patch  if  the  value  of  u 
and  v  for  that  point  is  known.  Note  that  in  equation 
(12)  Che  two  subscripts  of  r,  ru,  and  rv  refer  to 
particular  comers  of  a  patch  as  shown  in  Figures  9 
and  10. 

As  was  stated  earlier  not  only  are  the  Cartesian 
Coordinates  r  •  (x,y,s)  required  for  any  point  of  a 
patch  but  also  the  (u,v)  coordinates  of  that  point 
are  required.  These  can  easily  be  found  from  equa¬ 
tions  (10)  and  (11)  if  we  know  which  surface  patch 
the  point  is  in  and  the  values  of  p and  v for  the 
patch. 

Once  the  boundary  points  have  been  chosen  and 
the  blending  equation  (12)  has  been  defined  the 
differential  equations  (6)  can  be  solved  by  some 


Fig.  10  (;,0  coordinate  system  or.  surface. 


numerical  scheme,  in  particular,  poir:  successive 
over  relaxation  (S0R)  was  used.  The  initial  guess 
for  the  solution  was  taker,  as  a  simple  linear  inter¬ 
polation  between  the  values  on  the  inner  and  oucer 
curves.  For  the  evaluation  of  the  terms  in  the  dif¬ 
ferencial  equadons,  i.e.  to  evaluate  the  coefficients 
a,b,c,  J2,  AjU,  and  fi2v  in  (6)  the  blending  equation 
which  already  has  been  defined  is  used. 

RESULTS  AND  CONCLUSIONS 

Results  are  shown  for  prolate  ellipsoids  of 
various  thicknesses  surrounded  by  either  a  sphere  or 
an  ellipsoid,  Figs.  (11)  and  (12).  Also,  considered 
as  an  inner  body  is  a  thin  circular  cross  section 
fuselage,  Fig.  (13),  Intersecting  spheres  or  ellip¬ 
soids,  Figs.  (14)  and  (15),  and  a  wing-body  combina¬ 
tion.  For  the  prolate  ellipsoids  and  the  fuselage 
results  with  and  without  coordinate  concentrat ion  are 
shown. 

The  equations  presented  for  generating  a  three- 
dimensional  coordinate  system  between  two  arbitrary 
shaped  bodies  (or  boundary  surfaces)  can  be  solved 
numerically  and  give  good  results  for  all  cases 
attempted  in  this  research.  Inner  body  shapes  which 
are  simple  (e.g.  ellipsoids,  fuselages,  or  Intersect¬ 
ing  spheres  and  ellipsoids)  can  easily  be  combined 
with  sn  outer  body  to  yieltf  s  very  smooth  coordinate 
system  between  them. 


-  4  For  the  wing-body  combination  these  same  equa¬ 
tions  very  easily  adapt  and  generate  a  coordinate 
system  as  shown  in  Figures  16a  and  b.  This  applica¬ 
tion  requires  that  the  coordinates  on  the  surface  of 
the  body  be  transformed  by  equations  (6)  to  a  system 
where  one  coordinate  line  is  the  trace  of  the  wing 
and  the  body. 

At  this  cine  the  most  imposing  problem  in  the 
solution  of  equations  (A)  is  the  cuorespondence  of 
the  points  on  the  inner  and  outer  boundary  surfaces. 
As  was  discussed  in  reference  to  intersecting  ellip¬ 
soids  (Figures  2  and  3)  this  correspondence  may 
dictate  whether  ornot  the  solution  converges.  The 
"gap"  near  the  wing  in  Fig.  16  is  due  to  the  method 
used  to  correspond  the  points  on  the  body  to  those 
on  the  outer  boundary  surface.  Although  conceptually 
it  is  not  difficult  to  imagine  various  correspon¬ 
dences  of  the  points  on  the  two  boundary  surfaces 
it  canbe  quite  difficult  to  generate  this  correspon¬ 
dence  through  a  numerical  algorithm 
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Fig.  11 


Inner  body  a  thick  prolate  ellipsoid  with 
major  axis  2  and  minor  axis  ^3  surrounded 


by  a  sphere  of  radius  ■ 
contours  for  a  section 
for  all  (E,r.)  or  ( 1 , J) 
section  ;  -  const  (I  » 
(J.K)  values.  In  both 
in  ,  K  •  1. 


i.  (a)  Coordinate 
C  *  const  (K  »  11) 
values;  (b)  for  a 
1)  for  all  (r,t)  or 
cases  no  contraction 
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Tig.  12  Cases  (a)  »-  i  (b)  of  Fig.  11,  with  contract- 
t ion  in  • ,  K  »  1.05. 
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Fig.  13  Inner  body  a  thin  prolate  ellipsoid  uith 

major  axis  2  and  minor  axis  0.13  surrounded 
by  an  ellipsoid  with  major  axis  3  and  minor 
axis  2.4.  (a)  Coordinate  contours  for  a 

section  £  ■  constant  (I  «  1)  for  all  (n ,i) 
i  or  (J,K)  values;  (b)  for  a  section  (  •  const 

(K  -  13)  for  all  (£,n)  or  (I.J)  values.  In 
both  cases  no  contraction  in  n.  K  -  1. 
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(b) 

Fig.  15  (a)  Coordinates  on  surface  of  two  intersect 

ing  ellipsoids,  (b)  Coordinate  contours 
between  interesting  ellipsoids  and  a 
sphere  for  a  section  C  *  const  (I  *  1) 
for  all  (n.r)  0-  (J,K)  values.  Nc  con¬ 
traction  in  t,  K  *  1. 


<b) 


Fig.  14  (a)  Coordinates  on  surface  of  two  inter¬ 

secting  spheres. 

(b)  Coordinate  contours  between  intersecting 
spheres  and  an  outer  sphere  for  a  sect  ion 
£  •  const  (I  •  1)  for  all  (n,;)  or  (J,K) 
values.  No  contraction  in  n,  K  ■  1. 


Fig.  16  Wing-body  combination  where  wing  is  gen¬ 
erated  by  extending  lines  from  center  of 
body  through  points  on  surface  of  body  to 
outer  shpere.  A  represents  surfate  coordi¬ 
nates  on  wing.  £  represents  surface  coordi¬ 
nates  on  body.  C  represents  three-dimen¬ 
sional  coordinates  between  wing-body  with  x, 
y ,4s  coordinates,  (a)  top  view  of  wing- 
body  combination,  (b)  view  from  the  front 
of  wing-body  combination. 
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